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Foreword 


The  upper  ocean  is  an  area  of  vital  importance  to  naval  operations.  Changes 
in  the  density  structure  of  this  region  due  to  atmospheric  forcing  can 
significantly  affect  the  performance  of  both  acoustic  and  nonacoustic  sen¬ 
sors.  Despite  considerable  observation  and  study,  mixing  processes  in  the 
upper  ocean  are  still  far  from  being  fully  understood.  As  a  result,  current 
models  of  these  processes  are  fairly  empirical  and  somewhat  varied  in  their 
assumptions.  The  testing  of  mixed-layer  models  and  the  evaluation  of  their 
ability  to  reproduce  the  observed  response  of  the  upper  ocean  to  atmospheric 
forcing  is  the  most  practical  way  to  identify  accurate  and  useful  modeling 
parameterizations,  and  is  an  important  step  toward  improving  the  Navy’s 
upper-ocean  prediction  capability. 


A.  C.  Esau,  Captain,  USN 
Commanding  Officer,  NORDA 


Executive  summary 


Several  models  of  the  upper  mixed  layer  of  the  ocean  (Mellor-Yamada  Level 
2  and  2Vi ,  Niiler,  Garwood,  Price,  and  Therry-Lacarrere)  were  compared 
using  (a)  idealized  forcing  that  consisted  of  cases  of  wind  deepening,  heating, 
and  cooling;  (b)  data  from  Ocean  Stations  November  and  Papa;  (c)  data  taken 
during  the  Mixed-Layer  Experiment  (MILE);  and  (d)  data  taken  from  R/P 
FLIP  in  the  spring  of  1980  about  400  km  off  California.  Comparisons  with 
both  idealized  and  observed  forcing  show  the  differences  among  the  models 
to  be  significant.  Differences  are  especially  noticeable  for  the  deepening  of 
the  mixed  layer  in  the  fall  and  winter  due  to  wind  mixing  and  convection, 
and  for  the  shallowing  of  the  mixed  layer  during  light  winds  and  strong  heating. 
Although  evaluation  of  the  models  is  complicated  by  uncertainties  (primari¬ 
ly  with  regard  to  advective  effects  and  forcing),  the  results  suggest  certain 
deficiencies  in  some  of  the  mixing  parameterizations. 
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Testing  and  comparison  of  several  mixed-layer  models 


I.  Introduction 

Several  mixed -layer  (ML)  models— the  Mellor-Yamada 
Level  2  (MYL2)  and  Level  2Vi  (MYL2.5),  Niiler,  Gar¬ 
wood,  Price,  and  Therry-Lacarrere  (TL)  models — were 
tested  and  compared  for  their  suitability  for  use  in  opera¬ 
tional  ocean  forecasting.  These  models  are  fairly  represen¬ 
tative  of  the  current  range  of  oceanic  mixed-layer  models, 
and  are  among  the  most  commonly  used.  The  models  were 
evaluated  primarily  with  respect  to  their  accuracy  in 
reproducing  the  observed  response  of  the  mixed  layer.  The 
computational  requirements  of  the  models  on  a  VAX 
11/785,  the  computer  on  which  the  simulations  were  per¬ 
formed,  is  also  discussed. 

The  mixed-layer  models  described  in  the  scientific 
literature  vary  considerably  in  terms  of  their  formulation, 
complexity,  and  numerical  efficiency  (for  a  review  see  Gar¬ 
wood,  1979).  One  of  the  reasons  for  the  wide  variety  of 
mixed-layer  models  in  current  use  is  that  there  is  still  a 
fair  degree  of  uncertainty  regarding  the  primary 
mechanisms  responsible  for  turbulent  mixing  in  the  up¬ 
per  ocean.  Mixing  can  be  divided  into  two  broad  categories: 
forced  convection,  which  is  driven  by  the  wind;  and  free 
convection,  which  is  driven  by  vertical  density  instabilities. 
For  forced  convection,  the  possible  mechanisms  for  mix¬ 
ing  include  surface  and  internal  wave  effects,  shear  in¬ 
stabilities,  secondary  circulations  in  the  mixed  layer  (such 
as  Langmuir  cells),  and  interactions  among  these  processes. 
For  pure  free  convection,  mixing  is  caused  by  buoyancy 
instabilities.  However,  there  is  little  concurrence  as  to  how 
much  penetration  (erosion  of  the  stable  region  at  the  base 
of  the  mixed  layer  due  to  overshooting  of  the  buoyant 
plumes)  occurs  in  the  ocean.  In  addition,  since  some  wind 
is  usually  present,  free  convection  in  the  ocean  usually 
occurs  in  conjunction  with  some  degree  of  forced 
convection. 

Because  of  the  degree  of  uncertainty  regarding  the  validi¬ 
ty  of  the  parameterizations  employed  by  mixed-layer 
models,  probably  the  best  way  to  evaluate  them  for  use 
in  a  forecast  model  is  to  test  them  with  as  much  data  as 
possible.  Most  published  mixed-layer  studies  involve  the 
consideration  of  a  single  model,  usually  at  a  single  loca¬ 
tion  and  often  under  a  limited  range  of  conditions.  What 


we  have  attempted  here  is  to  intercompare  a  number  of 
models  with  data  from  several  locations  that  include  a  wide 
variety  of  conditions. 

Data  sets  suitable  for  testing  mixed-layer  models  are 
limited  in  number  because  they  must  include  reliable 
meteorological  observations  for  model  forcing  that  resolve 
diurnal  and  synoptic  variability,  as  well  as  adequate  sub¬ 
surface  data  for  model  initialization  and  verification.  Ad¬ 
ditionally,  it  is  best  if  advective  effects,  which  are  difficult 
to  estimate,  are  minimal. 

The  mixed-layer  models  were  tested  using  data  from 
weathership  stations  November  (N)  and  Papa  (P),  which 
are  located  in  a  relatively  quiescent  region  of  the  eastern 
North  Pacific;  data  taken  during  MILE  (the  Mixed-Layer 
Experiment),  conducted  near  Papa  in  the  fall  of  1977;  and 
data  taken  from  the  R/P  FLIP  about  400  km  off  Califor¬ 
nia  in  the  spring  of  1980,  which  include  high-resolution 
conductivity-temperature-depth  (CTD)  observations  of  the 
diurnal  mixed-layer  response.  The  different  data  sets  have 
their  particular  advantages.  The  ocean  station  data  allow 
testing  on  long  time  scales  under  a  wide  variety  of  seasonal¬ 
ly  varying  conditions.  The  data  from  the  mixed-layer  ex¬ 
periments  provide  more  accurate  and  extensive  observa¬ 
tions,  but  for  a  limited  period  of  time. 

In  addition  to  testing  with  data,  the  models  were  com¬ 
pared  for  some  simple,  idealized  forcing  conditions  that 
consisted  of  cases  of  wind-deepening,  heating,  and  cool 
ing.  Such  experiments  are  a  means  of  isolating  the  models’ 
performance  in  various  mixing  regimes.  They  provide  an 
indication  of  how  the  models  will  respond  to  real  forcing 
conditions  and  why  the  model  predictions  differ  from  each 
other  and  from  the  observed  ocean  response. 

II.  Description  of  mixed-layer 
models  tested 

Only  a  brief  description  of  the  models  is  provided  here. 
For  more  detail  the  reader  is  referred  to  the  papers  cited. 
Table  1  provides  a  summary  of  the  turbulence  mechanisms 
employed  by  the  models  for  mixing.  Note  that  although 
two  models  may  use  the  same  mechanism,  it  may  be  im¬ 
plemented  in  different  ways. 


Table  1.  Turbulence  mechanisms  used  by  the  mixed-layer  models.  Note  that  the  same  mechanism  may  be  implemented 
differently  in  different  models. 


model 

type* 

3 

surface  flux  (u;) 

shear 

convection 

penetrative  convection 

depth-dependent  dissipation 

MYL2 

diff 

no 

yes 

yes 

no 

yes 

MYL2.5 

diff 

not  effective 

yes 

yes 

very  little 

yes 

Niiler 

bulk 

yes 

only  when 
deepening 

yes 

yes 

not  for  u? 

Garwood 

bulk 

yes 

not  In  ver¬ 
sion  tested 

yes 

yes 

yes 

Price 

bulk 

no 

yes 

yes 

no 

yes 

TL 

diff 

not  effective 

yes 

yes 

very  little 

yes 

*  diff  *  differential  primitive  equation  mixed-layer  model, 
bulk  ■*  bulk  or  integrated  mixed-layer  model 


The  models  are  of  two  basic  types  based  on  their  method 
of  formulation,  differential  and  bulk.  The  MYL2,  MYL2.5, 
and  TL  models  are  differential  models  in  the  sense  that 
the  equations  for  the  conservation  of  momentum,  heat, 
salt,  and  turbulent  kinetic  energy  (TKE)  are  used  in  their 
primitive  form  and  are  not  integrated  over  the  mixed  layer. 
The  mixed  layer  for  these  models  is  defined  by  the  region 
where  the  local  TKE  is  large  enough  to  provide  a  certain 
minimum  level  of  vertical  mixing.  The  Niiler,  Garwood, 
and  Price  models  are  bulk  (or  integrated)  models  in  which 
the  mixed  layer  is  assumed  to  be  a  well-defined  layer  that 
is  uniform  in  temperature  and  salinity.  The  governing 
equations  for  the  bulk  models  are  obtained  by  integrating 
the  primitive  equations  over  the  depth  of  the  mixed  layer. 

The  MYL2  model  (Mellor  and  Yamada,  1974;  Mellor 
and  Durbin,  1975)  is  a  differential  model  in  which  the 
rate  of  vertical  mixing  is  described  by  eddy  coefficients 
that  depend  on  the  local  TKE.  The  TKE  is  calculated  us¬ 
ing  a  quasi-equilibrium  form  of  the  TKE  equation  in  which 
shear  production,  buoyancy  production,  and  dissipation 
fare  in  local  balance.  For  forced  convection,  the  MYL2 
(model  depends  only  on  shear  production  (the  instability 
!of  the  Ekman  current  shear)  to  generate  turbulence  for 
mixing.  The  depth  of  mixing  is  determined  by  the  depth 
where  the  shear  instability  is  stabilized  by  the  stratifica¬ 
tion,  which  is  where  where  the  Richardson  number  (the 
ratio  of  buoyancy  production  to  shear  production)  exceeds 
a  critical  value  (about  0.23  for  the  MYL2  model).  Since 
shear  production  is  the  only  source  of  TKE  where  the 
density  gradient  is  locally  stable,  the  MYL2  model  predicts 
no  convective  penetration,  only  a  convective  adjustment 
where  the  mixing  is  just  sufficient  to  relieve  the  density 
instability. 


The  MYL2.5  model  (Mellor  and  Yamada,  1977  and 
1982)  is  similar  to  the  MYL2  model,  except  that  the  TKE 
equation  includes  storage  and  diffusion  of  TKE.  The  dif¬ 
fusion  term  should  (ideally)  provide  a  means  for  supply¬ 
ing  TKE  to  the  base  of  the  mixed  layer  for  penetrative 
convection.  However,  the  gradient  approximation  that  is 
used  in  the  MYL2.5  model  for  this  term  has  been  generally 
criticized  as  being  inadequate  for  describing  this  process 
(Zeman  and  Lumley,  1976).  The  MYL2.5  model  also 
employs  slightly  different  turbulence  constants  that  result 
in  a  smaller  critical  Richardson  number  (about  0.20). 

The  Niiler  model  (Niiler,  1975;  Davis  et  al.,  1981b) 
is  a  bulk  mixed-layer  model.  The  depth  of  mixing  is 
governed  by  an  integrated  form  of  the  TKE  equation, 
which  includes  a  surface  flux  of  TKE  from  wave  effects 
(often  called  the  u*  term),  shear  production,  and  a  fixed 
amount  of  penetrative  convection.  When  the  mixed  layer 
is  deepening,  all  the  terms  in  the  TKE  equation  are  ac¬ 
counted  for.  However,  when  shallowing,  the  shear  term 
is  not  included  and  the  mixed-layer  depth  (MLD)  is  deter¬ 
mined  by  a  balance  between  the  surface  flux  of  TKE  and 
buoyancy  production. 

The  Garwood  model  (Garwood,  1977)  is  also  a  bulk 
mixed-layer  model.  This  model  includes  separate  budgets 
for  the  horizontal  and  vertical  components  of  TKE.  The 
TKE  equations  account  for  a  surface  flux  of  TKE  from 
waves  and  for  penetrative  convection.  The  model  includes 
a  dissipation  term  that  is  proportional  to  the  Coriolis 
parameter  to  account  for  the  effects  of  the  earth’s  rota 
tion  on  TKE  dissipation.  The  version  of  the  model  used 
here  does  not  include  shear  production  of  TKE  since  Gar¬ 
wood  (1977)  found  that  this  term  did  not  significantly  af¬ 
fect  the  depth  of  mixing  predicted  by  the  model.  Therefore, 


the  model  does  not  require  the  calculation  of  the  momen 
turn  budget,  which  reduces  the  amount  of  calculation  re¬ 
quired  for  a  one-dimensional  simulation. 

The  Price  model  (Price  et  al.,  1986)  depends  only  on 
shear  and  density  instability  for  mixing,  and  in  this  sense 
it  is  similar  to  the  MYL2  model  although  it  has  a  bulk 
formulation.  For  forced  convection,  the  depth  of  mixing 
is  governed  by  a  bulk  critical  Richardson  number  (0.65), 
and  for  free  convection  it  (like  the  MYL2  model)  provides 
for  a  convective  adjustment  but  no  convective  penetra¬ 
tion.  The  formulation  of  the  Price  model  differs  from  that 
of  the  other  bulk  models  in  that  the  temperature,  salini¬ 
ty,  and  momentum  profiles  at  the  base  of  the  mixed  layer 
are  smoothed  every  time  step  using  a  local  Richardson 
number  criteria.  This  has  the  effect  of  creating  a  smooth 
transition  region  at  the  base  of  the  mixed  layer  rather  than 
the  sharp  jump  or  discontinuity  common  to  most  other 
bulk  models. 

The  TL  model  (Therry  and  Lacarrere,  1983)  is  a  dif¬ 
ferential  model.  It  is  similar  to  the  MYL2.5  model,  but 
includes  a  parameterization  of  counter-gradient  buoyan¬ 
cy  fluxes  and  a  more  elaborate  formulation  of  the  diffu¬ 
sion  of  TKE  to  improve  the  description  of  free  convec¬ 
tion.  The  critical  Richardson  number  of  the  TL  model 
is  about  0.47. 

A  brief  word  is  in  order  about  the  tuning  or  adjust¬ 
ment  of  the  model  constants  for  the  models  tested  here. 
Since  optimum  values  for  many  of  the  constants  that  scale 
the  processes  parameterized  by  the  models  have  not  been 
firmly  established,  there  is  some  leeway  in  setting  the 
values  for  these  constants.  A  practical  goal  is  to  find  the 
values  for  the  model  constants  that  give  consistently  good 
predictions  in  different  regions  over  the  entire  range  of 
expected  conditions.  If  a  model  consistently  under-  or  over¬ 
predicts  the  mixed-layer  depth  by  a  similar  amount,  the 
model  constants  can  generally  be  adjusted  to  correct  the 
problem.  However,  having  to  use  different  model  constants 
to  get  good  results  in  different  regions  or  situations  sug 
gests  that  the  model  parameterizations  themselves  are 
lacking. 

The  constants  for  the  models  tested  here  were  defined 
as  reported  in  the  references  cited  with  the  following  ex¬ 
ceptions.  The  scaling  factor  for  the  turbulent  length  scale 
for  the  MYL2  and  MYL2.5  models  was  increased  from 
0.1  to  0.2  to  increase  the  rate  of  mixing  in  the  mixed 
layer.  The  values  for  the  Niiler  model  constants  were  as 
recommended  by  Davis  et  al.  (1981b)  from  simulations 
with  the  MILE  data.  The  constants  scaling  the  entrain¬ 
ment  flux  at  the  base  of  the  mixed  layer  and  the  “rota¬ 
tional”  dissipation  term  for  the  Garwood  model  were  set 


to  4.5  and  4.6,  respectively,  based  on  the  results  of  simula¬ 
tions  at  Stations  November  and  Papa. 

III.  Comparison  of  mixed-layer 
models  for  some  simple  forcing 
cases 

The  mixed-layer  models  were  tested  for  three  types  of 
idealized  forcing  to  compare  their  response  to  simple  forc¬ 
ing  conditions.  The  forcing  was  chosen  to  represent  three 
basic  regimes:  wind-deepening,  heating,  and  cooling. 

Wind-deepening  is  defined  to  occur  when  the  mixed 
layer  deepens  due  to  the  erosion  of  the  stably  stratified 
region  at  its  base  by  wind-generated  turbulence.  The  depth 
of  mixing  is  governed  by  the  strength  of  the  density 
stratification  and  the  magnitude  of  the  wind. 

In  the  heating  regime  the  mixed-layer  depth  is  governed 
by  a  balance  between  the  stabilizing  effect  of  surface  heating 
(or  a  positive  surface  buoyancy  flux)  and  the  effect  of  mix¬ 
ing  due  to  wind-generated  turbulence.  This  balance  governs 
the  mixed-layer  depth  during  periods  when  the  mixed  layer 
is  shallowing. 

In  the  cooling  regime,  a  net  surface  heat  loss  (or  a 
negative  surface  buoyancy  flux)  causes  the  mixed  layer 
to  deepen  due  to  convection.  Convection  usually  occurs 
in  the  mixed  layer  at  night  after  sunset,  and  is  the  domi¬ 
nant  mechanism  for  the  deepening  of  the  mixed  layer  in 
the  fall  and  winter  when  solar  heating  is  weak  and  there 
is  increased  evaporative  cooling  from  the  wind. 

The  density  stratification  and  forcing  for  these  tests  were 
chosen  to  be  representative  of  conditions  in  midlatitudes. 

A.  Wind-deepening 

For  the  wind-deepening  experiments  the  initial  condi¬ 
tions  were  taken  to  be  sea-surface  temperature  (SST)  equal 
to  24°C,  thermal  stratification  a  uniform  0.05°C/m,  salinity 
a  constant  35  ppt,  and  Ekman  velocity  zero.  The  latitude 
was  taken  to  be  29.91  °N  where  the  inertial  period  is  24 
hours.  For  each  of  the  models,  three  different  cases  of  wind- 
deepening  were  run  with  wind  stress  values  of  1,4,  and 
16  dynes/cm2.  The  surface  heat  and  salinity  fluxes  were 
zero. 

The  mixed-layer  depth  for  each  of  the  models  after  5 
days  is  listed  in  Table  2.  Figure  la  shows  mixed-layer  depth 
versus  time  for  the  4  dynes/cm2  wind-stress  case.  From 
Table  2  it  can  be  determined  that  the  predicted  mixed- 
layer  depth  is  approximately  proportional  to  the  square 
root  of  the  wind  stress  for  all  the  models. 
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Table  2.  Comparison  of  response  of  mixed-layer  models 
to  simple  forcing. 


A,  Wind-deepening:  Linear  initial  stratification  (0.05° C/m),  con¬ 
stant  wind  stress,  no  surface  heat  flux.  Table  values  are  mixed  - 
layer  depth  in  meters  after  5  days. 

Wind  stress  (dynes/cm2)  * 

1 

4 

16 

Mellor-Vamada  Level  2 

19 

41 

76 

Mellor-Yamada  Level  2V% 

18 

39 

76 

Niiler 

27 

54 

108 

Garwood 

26 

51 

104 

Price 

20 

44 

88 

Therry-Lacarrere 

22 

46 

87 

B.  Heating:  Small  initial  stratification,  constant  wind  stress  (1 
dyne/cm2),  constant  heat  flux.  Table  values  are  mixed-layer  depth 
in  meters  after  2  days 

Surface  heat  flux  (ly/day)  *= 

150 

600 

2400 

Mellor-Yamada  Level  2 

28 

16 

8 

Mellor-Yamada  Level  2  Vs 

28 

16 

8 

Niiler 

17 

4 

1 

Garwood 

33 

14 

4 

Price 

36 

18 

10 

Therry-Lacarrere 

38 

20 

10 

C.  Cooling:  Linear  initial  stratification  (0.05° C/m),  constant  wind 
stress  (1  dyne/cm2),  constant  heat  flux.  Table  values  are  mixed- 
layer  depth  in  meters  after  120  days. 

Surface  heat  flux  (ly/day)  =  - 

100 

—  200 

-300 

Mellor-Yamada  Level  2 

72 

105 

122 

Mellor-Yamada  Level  2 Vs 

72 

104 

122 

Niiler 

101 

127 

150 

Garwood 

90 

117 

140 

Price 

71 

102 

126 

Therry-Lacarrere 

76 

104 

126 

The  MYL2,  MYL2.5,  Price,  and  TL  models  show  a 
similar  response  to  the  wind-forced  deepening  of  the  mixed 
layer.  There  is  a  rapid  deepening  during  the  first  day,  and 
a  very  gradual  deepening  thereafter.  This  behavior  is 
characteristic  of  the  shear  instability  mechanism  that  these 
models  employ.  The  small,  inertial  period  oscillation  of 
the  mixed-layer  depth  for  these  models  is  caused  by  residual 
stratification  within  the  mixed  layer  and  the  fluctuating 
shear  at  the  base  of  the  mixed  layer  due  to  inertial  mo¬ 
tions.  The  primary  reason  for  the  difference  in  mixed-layer 
depth  between  the  MYL2,  MYL2.5,  and  TL  models  is 
due  to  the  difference  in  critical  Richardson  number— a 
higher  value  of  the  critical  Richardson  number  tends  to 
generate  a  deeper  mixed  layer.  The  diffusion  of  TKE  in 


the  MYL2.5  model  appears  to  have  little  effect  on  the 
mixed-layer  depth  for  these  cases  of  forced  convection. 

The  Niiler  model  responds  on  two  distinct  timescales. 
The  rapid  deepening  of  this  model  during  the  first  12  hours 
is  due  mainly  to  the  shear  production  term  in  the  TKE 
equation,  and  the  more  gradual  deepening  that  follows 
is  due  to  the  surface  TKE-flux  term  (Niiler,  1975).  It  is 
notable  that  the  Niiler  model  responds  much  more  quickly 
than  the  MYL2,  MYL2.5,  TL,  and  Price  models  in  spite 
of  the  fact  that  the  same  mechanism,  shear  production, 
is  primarily  responsible  for  the  deepening.  The  difference 
in  response  time  seems  to  be  due  to  the  instantaneous 
mixing  of  heat  and  momentum  within  the  mixed  layer 
that  is  inherent  in  the  bulk  formulation  of  the  Niiler  model. 
The  Garwood  model,  which  does  not  include  shear  pro¬ 
duction,  responds  a  little  less  quickly  than  the  Niiler  model 
during  the  first  12  hours,  but  then  begins  to  catch  up 
and  is  almost  as  deep  after  a  few  days.  The  persistent 
deepening  of  the  Niiler  and  Garwood  models,  relative  to 
the  others,  that  occurs  after  the  first  day  is  due  to  the 
surface  TKE  flux  term. 

B.  Heating 

For  the  heating  experiments,  the  initial  conditions  and 
model  parameters  were  as  for  the  wind-deepening  ex¬ 
periments  except  that  the  initial  temperature  from  the  sur¬ 
face  to  100  m  was  a  uniform  19°C.  The  models  were  forced 
with  a  constant  wind  stress  of  1  dyne/cm2  and  heat  fluxes 
of  150,  600,  and  2400  ly/day  (1  ly  =  1  cal/cm2).  (For 
reference,  the  mean  rate  of  heating  through  the  spring 
and  summer  in  midlatitudes  is  about  150  ly/day,  although 
daily  peaks  can  be  as  high  as  100  ly/hour.) 

The  mixed-layer  depths  predicted  for  the  heating  ex¬ 
periments  are  listed  in  Table  2.  The  Niiler  and  Garwood 
bulk  models  respond  instantly  to  a  change  in  the  surface 
forcing  when  shallowing  because  the  time  derivative  terms 
in  the  entrainment  equation  are  dropped  when  shallow¬ 
ing  and  the  new  mixed-layer  depth  is  solved  for  algebraical¬ 
ly.  The  response  of  the  other  models  to  surface  heating 
depends  on  the  strength  of  the  forcing,  but  is  on  the  order 
of  a  few  hours  as  increasing  stratification  suppresses  tur¬ 
bulent  mixing  below  a  certain  depth.  On  diurnal  and  longer 
timescales  this  difference  in  response  is  not  very  noticeable 
in  the  model  results. 

The  MYL2,  MYL2.5,  Price,  and  TL  models  predict 
similar  mixed-layer  depths.  As  in  the  wind-deepening  case, 
the  differences  in  mixed-layer  depth  are  primarily  governed 
by  the  differences  in  critical  Richardson  number,  which 
determines  the  balance  between  the  buoyancy  and  shear 
production  terms  in  the  TKE  equation.  It  should  be  pointed 


out  that  the  mixed  layer  for  the  differential  models  (MYL2, 
MYL2.5,  and  TL)  shows  a  significant  amount  of  stratifica¬ 
tion  when  surface  heating  is  present.  The  eddy  coefficients 
do  not  mix  the  mixed  layer  rapidly  enough  to  keep  the 
temperature  in  the  mixed  layer  uniform. 

For  these  heating  experiments  where  the  source  of  heat 
is  at  the  surface  (i.e.,  there  is  no  penetration  of  solar  radia¬ 
tion),  the  mixed-layer  depth  for  the  Niiler  model  is  pro¬ 
portional  to  the  Obukov  length  (the  ratio  of  the  surface 
flux  of  TKE  to  buoyancy  production).  The  Niiler  model 
shows  a  stronger  dependence  on  the  surface  buoyancy  flux 
than  the  other  models  where  the  mixed-layer  depth  for 
the  shallowing  case  depends  on  the  Ekman  length  scale, 
as  well  as  the  Obukov  length. 

The  Garwood  model  shows  a  stronger  response  to  sur¬ 
face  heating  than  the  MYL2,  MYL2.5,  Price,  and  TL 
models  when  the  surface  heat  flux  is  high,  but  does  not 
get  as  shallow  as  the  Niiler  model. 

An  idiosyncrasy  of  the  Price  model  is  that  the  mixed- 
layer  depth  undergoes  an  oscillation  when  the  mixed  layer 
is  shallowing  under  the  influence  of  a  constant  surface 
heat  flux  and  wind  stress.  The  effective  mixed-layer  depth 
is  determined  by  the  bottom  of  the  envelope  of  these  oscilla¬ 
tions,  since  this  is  the  region  that  is  kept  mixed.  The  oscilla¬ 
tion  is  inherent  in  the  algorithm  used  by  the  model. 
Although  a  bit  disturbing,  this  behavior  does  not  have 
an  evident  effect  on  the  predicted  density  and  momen¬ 
tum  profiles. 

C.  Cooling 

For  the  cooling  experiments  the  initial  conditions  and 
model  parameters  were  the  same  as  for  the  wind-deepening 
experiments.  The  models  were  forced  with  a  constant  wind 
stress  of  1  dyne/cm2  and  surface  heat  fluxes  of  - 100, 

—  200,  and  —  300  ly/day.  (The  mean  surface  heat  flux  in 
midlatitudes  from  October  through  January  is  on  the  order 
of  - 150  ly/day.)  Table  2  lists  the  mixed-layer  depths 
predicted  by  the  models  after  a  period  of  120  days,  and 
Figure  lb  shows  mixed-layer  depth  versus  time  for  the 

-  200  ly/day  case. 

Observations  of  convective  boundary  layers  in  the  at¬ 
mosphere  and  in  the  laboratory  indicate  that  convective 
penetration  is  usually  sufficient  to  generate  an  entrain¬ 
ment  buoyancy  flux  at  the  base  of  the  mixed  layer  that 
is  10-30%  of  (but  opposite  in  sign  to)  the  surface-driving 
buoyancy  flux  (Stull,  1976).  For  the  cooling  experiment 
of  Figure  lb  this  would  result  in  a  final  mixed-layer  depth 
of  108-125  m,  whereas  simple  convective  adjustment 
would  produce  a  final  mixed -layer  depth  of  99  m. 

The  shallowest  mixed-layer  depths  are  for  the.MYL2, 


MYL2.5,  Price,  and  TL  models,  which  give  similar  results. 
For  pure  convection  (no  wind  stress)  the  MYL2  and  Price 
models  can  provide  no  more  than  a  convective  adjustment 
of  the  density  profile.  Local  shear  production  provides  the 
only  source  of  TKE  for  mixing  for  these  models  when 
the  local  density  gradient  is  stable.  However,  for  a  mixed 
layer  deeper  than  70  or  80  m  the  shear  of  the  wind-driven 
current  is  too  weak  at  the  base  of  the  mixed  layer  to  pro¬ 
vide  much  entrainment  mixing  during  typical  wind 
conditions. 

For  the  MYL2.5  and  TL  models,  the  diffusion  term 
in  the  TKE  equation  can  supply  TKE  for  convective 
penetration.  However,  the  fact  that  the  mixed-layer  depths 
predicted  by  these  models  are  similar  to  those  predicted 
by  the  MYL2  and  Price  models  indicates  that  the 
parameterization  used  for  this  term  is  not  very  effective 
in  generating  convective  penetration.  This  was  noted  by 
Mellor  and  Yamada  (1977)  for  the  MYL2.5  model  in 
simulating  the  laboratory  convection  experiments  of  Willis 
and  Deardorff  (1974).  Therry  and  Lacarrere  (1983),  on 
the  other  hand,  obtained  a  negative  heat  flux  overshoot 
of  about  15%  with  their  model  in  simulations  of  the 
Wangara  atmospheric  boundary  layer  experiment.  The 
reason  for  this  discrepancy  with  the  results  here  is  not 
known. 

The  Niiler  model  is  adjusted  to  dissipate  83%  of  the 
TKE  generated  by  the  buoyancy  term  during  free  con¬ 
vection.  Hence,  the  TKE  used  for  entrainment  mixing 
is  17%  of  that  generated  by  the  buoyancy  term  plus  the 
surface  flux  of  TKE  whose  relative  importance  decreases 
with  increasing  mixed-layer  depth.  In  a  test  made  without 
the  surface  TKE  flux  term  with  a  surface  heat  flux  of 
—  200  ly/day,  the  final  mixed-layer  depth  was  reduced  from 
the  127  m  of  Figure  lb  to  118  m. 

For  free  convection  the  Garwood  model  dissipates  73% 
of  the  TKE  generated  by  the  buoyancy  term  when  the 
mixed  layer  is  very  shallow,  which  is  less  than  the  Niiler 
model.  However,  this  percentage  increases  as  the  mixed 
layer  deepens  (as  long  as  the  Coriolis  parameter  is  not 
zero),  and  in  Figure  lb  the  Garwood  model  deepens  less 
rapidly  than  the  Niiler  model  below  50  m. 

IV.  Mixed-layer  model  testing 
at  Ocean  Stations  November 
and  Papa 

Three-hourly  surface  meteorological  and  subsurface 
bathythermograph  (BT)  observations  have  been  collected 
at  the  ocean  weathership  stations  in  the  North  Atlantic 


and  North  Pacific  since  the  1940s.  This  is  probably  the 
best  available  data  for  testing  mixed-layer  models  on  synop¬ 
tic  to  interannual  time  scales.  Although  we  are  interested 
here  in  evaluating  mixed-layer  models  for  short  time-scale 
predictions,  there  are  a  number  of  advantages  in  testing 
the  models  on  longer  time  scales.  The  models  are  forced 
by  a  wide  range  of  seasonally  varying  conditions,  and  biases 
in  the  model  results  can  be  easier  to  discern  when  the 
simulations  are  carried  over  a  long  period  of  time. 

Stations  November  and  Papa  are  located  in  the  eastern 
North  Pacific  at  140°W,  30°N  and  145°W,  50°N,  respec¬ 
tively.  Data  for  the  years  1960-1970  were  investigated 
for  both  stations.  Because  of  large  gaps  in  the  BT  data 
at  Station  November  for  most  of  the  other  years,  1961 
was  chosen  for  testing  the  mixed-layer  models. 

A.  Initialization  and  forcing  from 
ocean  station  data 

The  BT  observations  were  used  to  initialize  the 
temperature  profiles  and  verify  the  model-predicted  SST 
and  mixed-layer  depth.  The  inclusion  of  salinity  in  the 
simulations  was  hindered  by  a  lack  of  salinity  and  precipita¬ 
tion  data.  At  Station  November  salinity  was  taken  to  be 
constant,  since  the  effect  of  salinity  on  the  density  stratifica¬ 
tion  of  the  upper  ocean  in  this  region  is  generally  small. 
At  Station  Papa  the  depth  of  mixing  in  winter  is  limited 
by  the  strong  halocline  which  exists  below  100  m  where 
the  thermal  stratification  becomes  relatively  weak.  There¬ 
fore,  to  prevent  an  unrealistically  deep  winter  mixed  layer, 
the  salinity  profile  at  Papa  was  initialized  from  climatology 
(Beatty,  1977).  Model  experiments  at  Station  Papa  in  which 
precipitation  was  estimated  showed  that  surface  salinity 
fluxes  generally  have  only  a  small  effect  on  vertical  mix¬ 
ing,  and  these  fluxes  were  set  to  zero  for  the  results  shown 
here. 

The  meteorological  data  taken  at  Stations  November 
and  Papa  include  SST,  air  and  wet  bulb  temperature,  wind 
speed  and  direction,  and  cloud  cover.  From  these  obser¬ 
vations,  standard  formulas  were  used  to  calculate  surface 
fluxes  of  momentum  and  heat  (see  Martin,  1985).  Solar 
radiation  was  calculated  using  the  Fritz  formula  for  clear  - 
sky  insolation  (List,  1958),  which  provides  a  proper  diur¬ 
nal  variation.  The  Reed  (1977)  doud  cover  correction  was 
used  for  Station  November  and  the  Tabata  (1964)  correc¬ 
tion  for  Papa.  The  albedo  was  taken  to  be  6%.  Seawater 
turbidity,  which  has  a  significant  effect  on  the  depth  of 
penetration  of  solar  radiation  and,  as  a  result,  the  develop¬ 
ment  of  the  upper  ocean  thermal  structure  in  spring  and 
summer  (Martin,  1985),  was  taken  to  be  Type  I  at 
November  and  Type  II  at  Papa  (Jerlov,  1976). 


B.  Results  of  mixed-layer  simulations  at 
Ocean  Stations  November  and  Papa  • 

The  observed  and  model-predicted  three-hourly  SST  and 
mixed-layer  depth  at  Stations  November  and  Papa  are 
shown  in  Figures  2  and  5,  and  Figures  3  and  6  show 
predicted  minus  observed  SST  for  the  two  stations.  The 
mixed-layer  depth  for  both  the  observations  and  the  simula¬ 
tions  is  defined  as  the  depth  at  which  the  temperature  ® 
becomes  0.1°C  less  than  the  SST.  A  different  criterion 
for  defining  mixed-layer  depth,  such  as  a  temperature  dif¬ 
ference  of  0.2°C,  misses  a  few  of  the  shallow,  short-term 
mixed  layers,  but  gives  a  similar  envelope  on  the  synop¬ 
tic  time  scale  for  both  the  observations  and  the  model 
results.  ® 

The  SST  determined  from  the  BT  data  at  November 
has  a  few  seemingly  spurious  features  that  should  be  kept 
in  mind  when  comparing  the  model  results  with  the  obser¬ 
vations.  The  more  notable  of  these  are  a  drop  of  about 
a  degree  during  the  middle  of  January,  a  sudden  increase 
of  about  a  degree  during  the  beginning  of  September  that  ® 
lasted  about  a  week  (this  feature  is  corrected  for  in  Figs. 

2  and  3),  and  three  abrupt  changes  near  the  beginning 
and  end  of  October.  These  SST  features  are  not  reflected 
in  the  mixed-layer  depth  or  surface  meteorological  obser¬ 
vations,  hence  they  are  not  caused  by  surface  forcing  and  ^ 
are  not  reproduced  by  the  models.  They  do  not  seem  to 
be  associated  with  a  change  of  the  reported  location  of 
the  observations,  but  they  do  frequently  coincide  with  a 
change  of  the  ship  on  station. 

Generally,  fluctuations  in  the  observed  SST  and  mixed - 
layer  depth  on  the  synoptic  time  scale  tend  to  be  well  £ 
correlated  with  the  surface  forcing  and,  as  a  result,  can 
be  seen  in  the  model  simulations.  In  fall  and  winter  the 
variability  predicted  by  the  models  is  less  than  observed. 
Especially  noticeable  are  the  shallow  mixed  layers  that  oc¬ 
cur  at  Station  November  during  the  winter  months.  These 
are  stabilized  by  temperature  jumps  of  up  to  0.8°C  and  £ 
appear  to  be  due  to  advection. 

All  the  models  predicted  a  premature,  overly  abrupt 
spring  shallowing  of  the  mixed  layer  at  Station  November. 

This  could  be  due  to  advection  or  to  an  overestimate  of 
the  surface  heat  flux  during  this  period.  The  time  of  the 
spring  shallowing  at  Papa,  however,  is  predicted  on  • 
schedule.  All  the  models  predicted  a  colder  SST  at  Papa 
than  observed  in  the  fall,  even  those  that  underpredicted 
the  mixed-layer  depth.  The  heat  content  of  the  mixed  layer 
at  Papa  at  this  time  cannot  be  accounted  for  by  the  sur¬ 
face  heat  flux,  and  must  have  been  due  to  a  net  advection 
of  heat  into  this  region.  ® 
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The  MYL2  and  MYL2.5  models  underpr edict  the 
mixed-layer  depth  at  both  November  and  Papa.  As  a  result 
the  predicted  SST  is  about  1.5°C  too  high  in  the  sum¬ 
mer.  Advective  effects  or  insufficient  mixing  are  the  most 
likely  causes  of  the  underpredicted  mixed-layer  depth.  It 
is  difficult  to  determine  if  advection  is  the  source  of  the 
discrepancy.  However,  evidence  that  the  MYL2  and 
MYL2.5  models  do  not  provide  enough  entrainment  mix¬ 
ing  is  that  the  temperature  gradient  at  the  base  of  the 
mixed  layer  is  often  significantly  less  than  observed, 
especially  in  late  summer  and  fall  at  Station  November 
(Fig.  4).  It  may  be  that  the  shear  production  parameterized 
by  these  models  is  insufficient  to  describe  forced  convec¬ 
tion  in  the  ocean  or  that  the  models  predict  too  little 
penetrative  convection. 

The  TL  model  gives  fairly  good  results  at  Papa,  and 
better  results  at  November  than  the  MYL2  and  MYL2.5 
models,  although  the  mixed-layer  depth  at  November  is 
still  shallower  than  observed.  The  improvement  over  the 
MYL2  and  MYL2.5  models  is  primarily  due  to  the  larger 
critical  Richardson  number  that  the  TL  model  employs. 

The  Niiler  model  does  fairly  well  in  simulating  the 
seasonal  evolution  of  SST  and  mixed-layer  depth  at  Sta¬ 
tion  November.  However,  the  predicted  mixed-layer  depth 
at  Papa  is  too  deep  and  the  summer  SST  is  underpredicted 
by  about  3°C.  Station  Papa  has  stronger  winds  than 
November  and  a  much  weaker  main  thermocline,  which 
provides  less  resistence  to  erosion  in  winter.  When  deepen¬ 
ing,  the  Niiler  model  takes  into  account  TKE  generated 
by  mean  shear,  the  surface  flux  from  waves,  and  convec¬ 
tion.  However,  when  shallowing,  only  the  surface  flux 
of  TKE  from  waves  balances  buoyancy  production.  Hence, 
the  scaling  of  the  surface  flux  must  be  sufficiently  large 
to  produce  realistic  mixed-layer  depths  when  shallowing, 
and  this  term  may  be  too  large  when  the  mixed  layer  is 
deepening. 

The  mixed-layer  depth  predicted  by  the  Garwood  model 
at  November  agrees  well  with  the  observations  through 
much  of  the  year,  but  is  a  little  shallower  than  observed 
in  the  summer,  particularly  during  periods  of  light  winds. 
At  Papa  the  agreement  between  the  Garwood  model  and 
the  observations  is  very  good  through  midsummer, 
although  the  predicted  mixed-layer  depth  is  slightly  deeper 
than  observed  in  the  fall.  The  agreement  was  not  improved 
by  additional  model  tuning;  e.g.,  decreasing  the  surface 
flux  of  TKE  or  increasing  the  dissipation  tends  to  make 
the  mixed  layer  too  shallow  in  summer.  However,  some 
of  the  discrepancy  may  be  due  to  the  net  advection  of 
heat  that  appears  to  have  occurred  at  Papa  in  .the  fall. 


The  Price  model  gives  results  similar  to  the  TL  model, 
a  fairly  good  prediction  of  mixed-layer  depth  at  Papa,  but 
too  shallow  at  November.  In  fact,  this  problem  was  found 
with  all  the  mixing  parameterizations  that  were  tried  which 
depend  on  shear-generated  turbulence  as  the  primary  mix¬ 
ing  mechanism  for  forced  convection.  It  may  be  that  shear 
is  insufficient  as  a  mixing  mechanism  for  forced  convec¬ 
tion,  or  that  some  other  error,  such  as  neglect  of  advec¬ 
tive  effects  at  November,  is  causing  the  discrepancy. 


V.  Mixed-layer  model  testing 
with  data  from  MILE 

MILE  was  conducted  near  ocean  station  Papa  (50°N, 
145°W)  in  the  eastern  North  Pacific  in  August  and 
September  1977.  This  experiment  was  specifically  con¬ 
ducted  to  investigate  the  upper  mixed  layer  and  to  pro¬ 
vide  data  against  which  models  could  be  tested.  Davis  et 
al.  (1981a, b)  analysed  some  of  the  data  from  MILE  and 
used  the  Niiler  model  to  simulate  the  response  of  the  mixed 
layer  during  the  experiment  to  local  atmospheric  forcing. 
They  noted  that  advective  contributions  to  the  upper  ocean 
heat  budget  during  MILE  were  generally  small,  but  not 
negligible,  and  parameterized  the  effects  of  advection  in 
their  model  simulations.  The  model  testing  with  the  MILE 
data  reported  here  is  based  on  their  work.  However,  no 
attempt  was  made  to  parameterize  advective  effects.  The 
local  heat  budget  in  the  model  simulations  is  due  only 
to  the  surface  heat  flux. 

A.  Model  initialization  and  forcing 
with  the  MILE  data 

At  MILE  mooring  1,  temperature  data  were  taken  at 
fixed  depths  from  5  to  175  m  at  intervals  of  112.5  sec 
from  August  19  to  September  7.  This  data  was  used  to 
initialize  and  validate  the  models. 

Meteorological  observations  were  taken  from  the  Cana¬ 
dian  weather  ship  Quadra  at  3-hour  intervals.  Standard  for¬ 
mulas  were  used  to  calculate  the  surface  wind  stress,  net 
longwave  radiation,  and  latent  and  sensible  heat  exchange 
from  the  wind  speed  and  direction,  air  temperature  and 
humidity,  and  cloud  cover.  Solar  radiation  was  measured 
directly  aboard  the  Quadra  at  hourly  intervals. 

B.  Results  of  mixed-layer  simulations 
with  the  MILE  data 

A  comparison  of  model-predicted  and  observed  5-m 
temperature  and  mixed-layer  depth  for  simulations  with 
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the  MILE  data  for  all  the  models  is  shown  in  Figure  8. 
The  temperature  at  5-m  depth  was  used  for  comparison 
because  there  are  no  reliable  measurements  of  the  actual 
SST  during  MILE.  The  mixed-layer  depth  is  again  defined 
as  the  depth  at  which  the  temperature  is  0.1°C  less  than 
the  surface  temperature. 

The  observed  5-m  temperature  during  MILE  shows  the 
steady  warming  of  the  ocean  that  was  still  occurring  late 
in  the  summer.  The  steady  warming  was  modulated  by 
a  diurnal  fluctuation  of  0.2-0. 3°C  caused  by  daytime  solar 
heating  and  nighttime  cooling.  A  cooling  event  occurred 
on  August  22  and  23  caused  by  an  increase  in  the  winds, 
and  a  smaller  cooling  event  occurred  at  the  end  of  August. 

The  models  generally  predict  the  trend  and  the  diurnal 
fluctuation  of  the  5-m  temperature  fairly  well.  None  of 
the  models  predicted  as  much  cooling  as  occurred  on 
August  22  and  23,  even  those  that  predicted  the  observed 
amount  of  mixed-layer  deepening,  which  indicates  that 
some  of  the  cooling  at  that  time  was  due  to  advection. 

The  5-m  temperature  and  mixed-layer  depth  predicted 
by  the  MYL2  model  agree  fairly  well  with  the  observa¬ 
tions.  The  major  discrepancy  is  that  the  model  under¬ 
estimates  the  amount  of  cooling  that  occurred  on  August 
22  and  23.  Part  of  the  discrepancy  is  due  to  an  under¬ 
estimate  of  the  mixed-layer  deepening,  and  part  to  advec¬ 
tion  not  accounted  for  in  the  model  simulations.  As  a  result 
of  underpredicting  the  cooling  on  August  22  and  23,  the 
5-m  temperature  predicted  by  the  MYL2  model  remains 
biased  by  about  0.4°C  for  the  rest  of  the  experiment.  For 
a  simulation  initialized  on  August  25,  after  the  cooling 
event,  the  MYL2  model  tracks  the  observed  5-m 
temperature  very  well  (Fig.  9).  The  MYL2.5  model,  as 
at  Stations  November  and  Papa,  gives  results  similar  to 
the  MYL2  model,  but  with  a  slightly  shallower  and,  as 
a  result,  warmer  mixed  layer. 

The  Niiler  model  gives  a  good  prediction  of  mixed-layer 
depth  during  the  cooling  event  of  August  22  and  23,  but 
afterward  generally  overpredicts  the  mixed-layer  depth  and 
the  mixed  layer  warms  less  quickly  than  observed.  The 
Garwood,  Price,  and  TL  models  give  similar  results,  which 
lie  between  those  of  the  models  discussed  previously,  and 
show  fairly  good  overall  agreement  with  the  observed  5-m 
temperature  and  mixed-layer  depth. 

VI.  Mixed-layer  model  testing 
with  data  from  R/P  FLIP 

A  useful  data  set  for  mixed-layer  model  testing  was  ob¬ 
tained  from  R/P  FLIP  in  April  and  May  1980  about  400 


km  west  of  San  Diego,  California,  at  about  30°N.  The 
frequent  CTD  observations  that  were  taken  provide  a 
detailed  look  at  the  response  of  the  mixed  layer  to  diur¬ 
nal  solar  forcing  under  a  variety  of  wind  conditions.  Price 
et  al.  (1986)  reported  on  the  observations  and  on  simula¬ 
tions  of  the  mixed-layer  response  using  the  Price  mixed- 
layer  model.  The  mixed-layer  simulations  performed  for 
model  testing  here  were  based  on  their  study. 

A.  Model  initialization  and  forcing 
with  FLIP  data 

Temperatures  taken  by  a  profiling  CTD  every  120  sec 
with  a  vertical  resolution  of  1  m  were  used  to  initialize 
and  validate  the  models.  The  uppermost  temperature 
measured  was  at  a  depth  of  2  m. 

Surface  wind  stress,  net  longwave  radiation,  and  latent 
and  sensible  heat  exchange  were  calculated  from  half-hourly 
meteorological  observations  using  standard  formulas.  Solar 
radiation  was  measured  directly  with  an  Eppley 
pyranometer. 

There  were  significant  advective  effects  during  most  of 
the  experiment.  However,  for  a  4-day  period  between  May 
7  and  11  (Julian  day  128-132),  advective  effects  were 
minimal  (Price  et  al.,  1986),  and  this  period  was  used  for 
model  testing. 

B.  Results  of  mixed-layer  simulations 
with  FLIP  data 

Figure  10  shows  hourly  observed  and  model-predicted 
temperature  profiles  for  May  9  and  10.  During  the 
preceding  two  days  the  winds  were  moderate  (8-12  m/sec) 
and  the  mixed  layer  stayed  fairly  well  mixed  down  to  about 
40  m  near  the  top  of  the  seasonal  thermocline.  During 
the  two  days  shown  the  winds  were  light  (3-7  m/sec)  and 
a  shallow  mixed  layer  formed  in  response  to  diurnal  solar 
heating. 

Notable  features  in  the  observed  profiles  are  that  (1) 
there  is  almost  always  a  fairly  well-mixed  surface  layer, 
i.e.,  the  temperature  in  the  surface  layer  is  quite  uniform; 
and  (2)  there  is  a  fairly  smooth  transition  region  at  the 
base  of  the  mixed  layer.  The  uniformity  of  the  mixed  layer 
may  be  due  to  the  Langmuir-type  circulations  observed 
(Weller  et  al.,  1985),  since  such  overturning  circulations 
can  mix  very  rapidly.  The  transition  region  at  the  base 
of  the  well-mixed  layer  appears  to  be  a  region  where  the 
mixing  rates  are  considerably  reduced,  but  are  still 
significant. 

The  MYL2  and  MYL2.5  models  again  show  similar 
results.  Consistent  with  the  findings  at  Stations  November 
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and  Papa,  the  mixed  layer  is  not  quite  as  deep  as  observed. 
Also  apparent  here,  when  compared  to  the  highly  resolved 
CTD  profiles,  is  that  there  is  often  more  stratification  in 
the  mixed  layer  than  observed,  which  suggests  that  the 
rate  of  mixing  in  the  mixed  layer  predicted  by  these  models 
is  too  slow.  This  occurs  both  during  the  day  when  there 
is  heating,  and  at  night  when  there  is  cooling  and  con¬ 
vection.  These  turbulent  diffusion  models  were  developed 
to  describe  turbulent  boundary  layers  over  rigid  surfaces, 
and  do  not  account  for  circulations  that  might  have  the 
same  scale  as  the  boundary  layer  and  cause  more  rapid 
mixing  over  the  mixed  layer  than  smaller  scale  turbulent 
eddies.  The  diffusion  models  do,  however,  provide  the  kind 
of  smooth  transition  region  at  the  base  of  the  mixed  layer 
seen  in  the  FLIP  CTD  observations. 

The  Niiler  model  predicts  a  shallower  mixed  layer  than 
observed  when  there  is  strong  heating  and  light  winds. 
From  the  idealized  forcing  experiments  we  know  that  this 
model  responds  much  more  strongly  than  the  others  to 
a  high  surface  heat  flux.  ‘The  FLIP  data  suggest  that  this 
response  is  too  strong.  The  Garwood  model  responds  less 
strongly  than  the  Niiler  model  to  these  conditions,  but 
still  more  strongly  than  observed.  Also,  the  bulk 
parameterization  used  in  the  Niiler  and  Garwood  models 
gives  a  sharper  gradient  at  the  base  of  the  mixed  layer 
than  observed  in  the  FLIP  CTD  profiles. 

The  Price  model  gives  a  good  simulation  of  the  FLIP 
data.  The  SST  and  mixed-layer  depth  agree  fairly  well  with 
the  observations,  and  the  bulk  Richardson  number  deepen¬ 
ing  criteria,  combined  with  a  local  Richardson  number 
criteria  for  smoothing  the  profile  at  the  base  of  the  mixed 
layer,  generates  a  realistic-looking  (compared  to  the  FLIP 
CTD  observations)  temperature  profile. 

The  mixed  layer  predicted  by  the  TL  model  is  slightly 
deeper  than  predicted  with  the  MYL2  and  MYL2.5 
models,  which  seems  to  give  better  agreement  with  the 
observations.  But,  as  with  the  other  differential  models, 
the  mixed  layer  is  not  as  well  mixed  as  observed. 

VTL  Computational  requirements 
for  the  mixed-layer  models 

For  a  mixed-layer  model  applied  at  a  single  location  for 
a  forecast  of  short  duration,  the  computer  time  and  storage 
requirements  of  the  models  considered  here  are  relative¬ 
ly  small.  However,  if  a  mixed-layer  model  is  run  for  a 
long  period  of  time  or  is  applied  at  a  large  number  of  points, 
the  time  and/or  storage  requirements  can  become 
significant. 


Table  3  lists  the  central  processor  (cp)  time  required 
for  a  1-day  forecast  at  a  single  location  with  the  various 
models.  The  timings  were  done  on  a  VAX  11/785,  a  scalar 
(nonvectorizing)  computer  with  a  calculation  speed  of  about 
300,000  operations  per  second.  The  performance  of  the 
differential  models  (MYL2,  MYL2.5,  TL)  relative  to  the 
bulk  models  would  be  considerably  better  on  a  vectoriz¬ 
ing  computer  because  of  the  high  degree  of  vectorizabili- 
ty  of  the  differential  models. 


Table  3.  Central  processor  (cp)  time  used  by  mixed-layer 
models  for  1-day  simulation  on  VAX  11/785  computer. 


model 

timestep 

cp  time 

Mellor-Yamada  Level  2 

10  min 

3.0  sec 

Mellor-Yamada  Level  2% 

10  min 

7.0  sec 

Niiler 

1  hr 

0.4  sec 

Garwood 

1  hr 

0.1  sec 

Price 

1  hr 

1 .0  sec 

Therry-Lacarrere 

10  min 

5.0  sec 

The  cp  times  in  Table  3  must  be  qualified  by  the  follow¬ 
ing  comments.  First,  the  numerical  algorithms  used  were 
not  necessarily  developed  for  maximum  efficiency.  The 
differential  models  (MYL2,  MYL2.5,  and  TL)  and  the 
Niiler  bulk  model  were  programmed  inhouse  at  NORDA 
for  research  purposes  with  the  emphasis  on  ease  of  use 
and  modification,  rather  than  on  maximum  speed.  The 
algorithms  for  the  Garwood  and  Price  bulk  models  were 
obtained  from  their  authors.  The  differential  models  have 
the  potential  to  be  speeded  up  by  a  factor  of  2  or  3  by 
using  more  efficient  algorithms.  The  Niiler  and  Price  bulk 
models  have  some  potential  for  being  speeded  up  as  well. 
The  Garwood  model,  however,  appears  to  be  fairly  effi¬ 
cient  in  its  present  form. 

Second,  the  computer  time  required  for  a  model  simula¬ 
tion  depends  on  grid  resolution,  timestep,  number  of  grid 
points  and,  to  some  degree,  the  particular  circumstances 
of  the  simulation.  The  times  listed  in  Table  3  are  based 
on  an  average  of  the  times  required  for  the  simulations 
reported  here  with  the  November,  Papa,  MILE,  and  FLIP 
data,  and  thus  represent  a  wide  range  of  conditions.  The 
times  are  also  based  on  the  use  of  a  stretched  grid  of  45 
points  from  the  surface  to  300-m  depth  with  a  resolution 
of  2  m  at  the  surface.  Such  a  grid  has  adequate  resolution 
and  depth  for  most  applications.  The  model  timestep  used 
in  determining  the  model  cp  times  was  10  min  for  the 
differential  models  and  1  hr  for  the  bulk  models.  The 
timestep  for  the  bulk  models  is  limited  by  the  need  to 
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resolve  the  variability  in  the  surface  forcing  and  to  resolve 
the  bulk  or  integrated  response  of  the  mixed  layer  to  that 
forcing.  The  timestep  for  the  differential  models,  on  the 
other  hand,  is  further  limited  by  the  need  to  provide  an 
accurate  treatment  of  the  vertical  diffusion  terms.  Although 
the  diffusion  terms  are  treated  implicitly,  instabilities  can 
develop  if  the  timestep  is  too  much  larger  than  the  dif¬ 
fusive  timescale  based  on  the  grid  spacing  (the  diffusive 
timescale  equals  the  square  of  the  grid  spacing  divided 
by  the  eddy  diffusion  coefficient). 

The  cp  times  listed  in  Table  3  range  from  7  to  0.1  sec 
per  model  day,  with  the  MYL2.5  model  being  the  slowest 
and  the  Garwood  model  the  fastest.  The  differential  models 
generally  require  more  computation  than  the  bulk  models 
because  of  the  need  to  use  a  smaller  timestep.  The  dif¬ 
ferences  in  speed  among  the  differential  models  is  due  to 
the  increased  computation  required  in  the  eddy  coefficient 
algorithms  for  the  TL  and  MYL2.5  models  where  the 
mean  TKE  equation  must  be  solved.  Of  the  bulk  models, 
the  Price  model  is  somewhat  slow  because  of  the  time 
consumed  by  the  relaxation  scheme  used  to  adjust  the 
transition  region  at  the  base  of  the  mixed  layer.  The  speed 
of  the  Garwood  model  is  due  to  its  simplicity  (the  momen¬ 
tum  budget  is  not  required)  and  the  use  of  efficient 
algorithms  to  calculate  changes  in  mixed-layer  depth. 

The  computational  disadvantage  of  differential  models 
with  respect  to  bulk  models  (due  to  the  smaller  timestep 
required)  is  somewhat  offset  on  large  “vectorizing”  com¬ 
puters  because  of  the  high  degree  of  vectorizability  of  the 
differential  models;  i.e.,  the  same  computations  are  per¬ 
formed  at  each  grid  point  each  time  step.  The  algorithms 
for  bulk  models  tend  to  be  difficult  to  vectorize,  since  dif¬ 
ferent  sequences  of  calculations  are  performed  at  different 
points. 

For  applications  requiring  more  (or  less)  resolution, 
changing  the  grid  resolution  affects  the  time  requirements 
of  the  models  differently.  For  the  bulk  models,  increasing 
the  grid  resolution  results  in  a  somewhat  less  than  pro¬ 
portional  increase  in  the  cp  time  required;  e.g.,  doubling 
the  grid  resolution  and  the  number  of  grid  points  may 
increase  the  cp  time  from  10  to  100%  depending  on  the 
amount  of  grid-dependent  calculation  involved.  For  the 
differential  models,  a  doubling  of  the  grid  resolution 
generally  requires  a  four-fold  decrease  in  the  time  step 
to  maintain  the  accuracy  of  the  vertical  diffusion  terms. 
Since  decreasing  the  timestep  increases  the  cp  time  pro¬ 
portionally  (true  for  all  the  models  tested  here),  doubling 
the  grid  resolution  for  the  differential  models  results  in 
about  a  6-  to  8-fold  increase  in  the  cp  time  required. 


VIII.  Summary 

Table  4  summarizes  the  results  from  the  mixed-layer 
model  testing  with  the  November,  Papa,  MILE,  and  FLIP 
data,  and  Table  5  summarizes  the  general  pros  and  cons 
of  the  mixed-layer  models. 


Table  4.  Specific  model  results  from  mixed-layer  model 
testing. 


model 

results 

MYL2 

Underpredicts  MLD  at  both  N  and  P  ML  not  as  well- 
mixed  or  as  deep  as  observed  with  FLIP  data. 

MYL2  5 

Gives  results  similar  to  MYL2.  Gradient  diffusion  of 
TKE  provides  little  modification  of  predicted  MLD, 
even  for  the  convective  case. 

Niiler 

Good  prediction  of  MLD  at  N.  Overpredicts  MLD  at 

P  where  main  thermodine  is  weak,  due  to  insufficient 
dissipation  when  ML  is  deep.  Overly  strong  response 
to  high  heat  flux/light  winds  with  FLIP  data  Gradient 
at  the  base  of  the  ML  generally  sharper  than 
observed. 

Garwood 

Fairly  good  predictions  at  N  and  P,  although  model 
tends  to  shallow  too  much  in  midsummer,  especial¬ 
ly  when  winds  are  light.  Response  to  light  winds/high 
heat  flux  with  FLIP  data  may  be  too  strong.  Gradient 
at  base  of  mixed  layer  generally  sharper  than 
observed. 

Price 

MLD  not  as  deep  as  observed  at  N.  Good  predic¬ 
tion  at  P  Good  simulation  of  FUP  data. 

TL 

MLD  too  shallow  at  N,  fairly  good  at  P.  ML  not  as  well- 
mixed  as  observed  with  FLIP  data 

The  MYL2  and  MYL2.5  models  tend  to  underpredict 
the  mixed-layer  depth.  The  TL  model,  which  has  a  larger 
critical  Richardson  number,  predicts  a  slightly  deeper 
mixed  layer  that  generally  agrees  better  with  the  obser¬ 
vations  (although  the  mixed-layer  depth  is  still  under- 
predicted  at  November).  With  a  change  in  turbulence  con¬ 
stants  to  increase  the  critical  Richardson  number,  the 
MYL2  and  MYL2.5  models  would  also  give  better  agree¬ 
ment  with  the  observed  mixed-layer  depth. 

The  mixed  layer  predicted  with  all  the  differential  models 
(MYL2,  MYL2.5,  and  TL)  was  not  as  well  mixed  (uniform) 
as  shown  by  the  high-resolution  CTD  profiles  taken  from 
FLIP.  This  evidence  suggests  that  the  eddy  coefficients 
predicted  for  the  mixed  layer  by  these  models  are  not  giv¬ 
ing  as  rapid  mixing  as  actually  occurs.  This  deficiency  could 
be  partially  alleviated  by  increasing  the  turbulent  mixing 
scale  in  these  models  (the  resulting  increase  of  the  eddy 
coefficients  would  require  a  decrease  in  the  timestep  to 
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Table  5.  Pros  and  cons  of  mixed-layer  models. 


model 

pro 

con 

MYL2 

Smooth  transition  at  ML  base 

Potential  to  predict  current  shears. 

Very  vectorizabla 

ML  too  shallow. 

ML  not  sufficiently  well-mixed. 

Not  very  fast. 

MYL2.5 

Same  as  MYL2. 

TKE  diffusion  helps  stabilize  numerics. 

Same  as  MYL2. 

Slowest  of  models  tested 

Niiler 

Large  number  of  TKE  generation 
mechanisms  provides  flexibility. 

Insufficient  dissipation  for  deep  ML 

Shear  not  accounted  for  when  shallowing. 

ML  too  shallow  with  high  heating/Hght  winds. 

Discontinuity  at  ML  base 

Garwood 

Fast. 

No  momentum  budget  needed. 

Shear  mixing  not  explicitly  accounted  for. 

ML  may  be  too  shallow  with  high  heating/light  winds. 
Discontinuity  at  ML  base. 

Price 

Simple. 

Smooth  transition  at  ML  base 

Constant  positive  heat  flux  causes  oscillation  of  MLD. 
Somewhat  slow  for  bulk  ML  model. 

May  not  mix  deep  enough  in  some  situations. 

TL 

Same  as  MYL2.5 

ML  not  sufficiently  wen-mixed. 

May  not  mix  deep  enough  in  some  situations. 

maintain  the  accuracy  of  the  diffusion  terms).  The  dif¬ 
ferential  models  do  predict  one  of  the  significant  features 
of  the  FLIP  temperature  profiles— a  smooth  transition  (gra¬ 
dient)  at  the  base  of  the  mixed  layer. 

The  Niiler  bulk  model  gives  a  good  simulation  of  SST 
and  mixed-layer  depth  at  November,  but  predicts  too  deep 
a  mixed  layer  at  Papa  and  with  the  MILE  data.  The  FLIP 
data,  in  turn,  suggest  that  this  model  predicts  too  shallow 
a  mixed  layer  when  the  winds  are  light  and  the  heating 
is  strong.  These  problems  can  be  traced  to  several  short¬ 
comings.  There  is  not  enough  dissipation  of  the  surface 
flux  of  TKE  when  the  mixed  layer  is  deep;  and  the  pro¬ 
duction  of  TKE  from  shear  instability,  which  is  taken  into 
account  when  the  mixed  layer  is  deepening,  is  not  ac¬ 
counted  for  when  the  mixed  layer  is  shallowing.  The  FLIP 
observations  also  indicate  that  the  gradient  predicted  by 
the  Niiler  model  at  the  base  of  the  mixed  layer  is  overly 
sharp,  a  characteristic  of  most  bulk  mixed-layer  models. 

The  Garwood  bulk  model  gave  the  best  overall  agree¬ 
ment  for  the  annual  simulations  at  both  November  and 
Papa.  The  model  did  benefit  from  being  tuned  somewhat 
to  these  two  data  sets,  although  some  improvement  could 
probably  be  achieved  with  additional  experimentation.  At 
November,  the  model  was  a  little  too  shallow  in  the  sum¬ 
mer,  and  at  Papa  the  model  tended  to  be  too  shallow  in 
the  summer  or  too  deep  in  the  fall  depending  on  model 
tuning.  The  FLIP  data  indicate  that  the  model  may  shallow 
too  much  when  strong  heating  and  light  winds  occur, 
though  not  as  severely  as  the  Niiler  model.  As  with  the 


Niiler  bulk  model,  the  Garwood  model  predicts  a  sharper 
gradient  at  the  base  of  the  mixed  layer  than  is  observed. 

The  Price  model  does  not  predict  as  deep  a  mixed  layer 
as  observed  at  November,  but  gives  fairly  good  predic¬ 
tions  with  the  other  data  sets.  It  gives  the  best  simulation 
of  the  CTD  observations  of  the  diurnal  mixed-layer  cycle 
taken  from  FLIP,  where  the  combination  of  using  a  bulk 
Richardson  number  to  deepen  the  mixed  layer  combined 
with  a  gradient  Richardson  number  to  diffuse  the  inter¬ 
face  generates  temperature  profiles  that  have  much  of  the 
character  of  those  observed.  The  relaxation  scheme  used 
for  diffusing  the  interface,  however,  results  in  the  Price 
model  being  relatively  slow  for  a  bulk  model.  When  forced 
with  a  steady,  positive  surface  heat  flux  and  wind  stress 
the  Price  model’s  mixed-layer  depth  undergoes  an  oscilla¬ 
tion  that  is  inherent  in  the  algorithm  used.  However,  the 
oscillation  does  not  appear  to  significantly  affect  the 
predicted  profiles. 

The  running  times  per  model  day  for  the  models  (cen¬ 
tral  processor  time  on  a  VAX  11/785)  ranged  from  7  sec 
for  the  MYL2.5  model  to  0.1  sec  for  the  Garwood  model. 
It  must  be  noted  that  the  models  were  not  necessarily 
optimized  for  maximum  speed.  However,  differential  (dif¬ 
fusion)  models  generally  require  more  computation  than 
bulk  models  because  of  the  need  to  use  a  smaller  timestep 
to  maintain  the  accuracy  and  stability  of  the  vertical  dif¬ 
fusion  terms.  This  disadvantage  would  be  somewhat  off¬ 
set  on  a  vectorizing  computer  because  of  the  high  degree 
of  vectorizability  of  the  differential  models. 
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Figure  1.  Mixed-layer  depth  versus  time  for  (a)  wind-deepening  experiment ,  and  (b)  cooling  experiment  with  MYL2  (solid 
line),  MYL2.5  (dotted  line),  Niiler  (dashed  line),  Garwood  (dash-dot  line),  Price  (dash-dot-dot  line),  and  TL  (dash-dot-dot-dot 
line)  models. 
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Figure  2a.  Observed  and  predicted  SST  and  MLD  at  Ocean  Station  November  for  the  year  196L  Predicted  results  are  for 
the  MYL2,  MYL2.5 ,  and  Niiler  models . 
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Figure  2b  SST  and  MLD  for  Ocean  Station  November  simulations  with  the  Garwood \  Price ,  and  TL  models. 
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Figure  3a.  Difference  between  predicted  and  observed  SST  at  November  with  the  MYL2,  MYL2.5,  and  Niiler  models 
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Figure  3b.  Difference  between  predicted  and  observed  SST  at  November  with  the  Garwood,  Price,  and  TL  models. 
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Figure  4a.  Observed  and  predicted  profiles  at  November  for  a  4  -day  period  near  the  end  of  October  (October  21-24).  Predicted 
profiles  are  for  the  MYL2,  MYL2.5,  and  Niiler  models.  The  horizontal  scale  is  1  day  =  10°C. 
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Figure  4b  Observed  and  predicted  profiles  at  November  for  a  4 -day  period  near  the  end  of  October  (October  21-24).  Predicted 
profiles  are  for  the  Garwood ,  Price ,  and  TL  models .  The  horizontal  scale  is  1  day  =  10°C. 
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{Figure  5a.  Observed  and  predicted  SST  and  MLD  at  Ocean  Station  Papa  for  the  vear  1961 
MYL2,  MYL2.5,  and  Niiler  models. 
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Figure  5b.  SST  and  MLD  for  Ocean  Station  Papa  simulations  with  the  Garwood,  Price,  and  TL  models. 
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Figure  6a.  Difference  between  predicted  and  observed  SST  at  Papa  with  the  MYL2,  MYL2.5 ,  and  Niiler  models. 
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Figure  6b.  Difference  between  predicted  and  observed  SST  at  Papa  with  the  Garwood,  Price,  and  TL  models. 
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Figure  la.  Observed  and  predicted  profiles  at  Papa  for  a  5-day  period  from  August  29  to  September  2.  Predicted  profiles  are 
for  the  MYL2,  MYL2.5,  and  Niiler  models.  The  horizontal  scale  is  1  day  =  10°C. 
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Figure  7b .  Observed  and  predicted  profiles  at  Papa  for  a  5 -day  period  from  August  29  to  September  2.  Predicted  profiles 
are  for  the  Garwood ,  Price ,  and  TL  models .  Tfte  horizontal  scale  is  1  day  =  10  °C. 
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Figure  8a.  SST  and  MLD  for  simulation  of  MILE  data  with  MYL2 ,  MYL2.5 ,  and  Niiler  models.  The  solid  line  shows  the 
model  results ;  the  dashed  line  shows  the  observations. 
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Figure  8b  SST  and  MLD  for  simulation  of  MILE  data  with  Garwood ,  Price ,  and  TL  models .  The  solid  line  shows  the  model 
results;  the  dashed  line  shows  the  observations. 


21 


M0DEL  AND  0BSERVED  5M  TEMPERATURE  F0R  MILE 


AUGUST  SEPTEMBER 


0BSERVED  AND  PREDICTED  MLD  F0R  MILE 


AUGUST  SEPTEMBER 


Figure  9  SST  and  MLD  for  simulation  of  MILE  data  with  MYL2  model  initialized  at  OOZ  August  25,  after  the  August 
22-23  storm.  The  solid  line  shows  the  model  results ;  the  dashed  line  shows  the  observations. 
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Figure  10a.  Observed  and  predicted  temperature  profiles  for  simulation  of  the  data  from  R/P  FLIP for  May  9-10,  1980.  Predicted 
profiles  are  for  the  MYL2,  MYL2.5,  and  Niiler  models.  The  horizontal  scale  is  1  day  =  1.43°C. 
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Figure  10b.  Observed  and  predicted  temperature  profiles  for  simulation  of  the  data  from  R/P  FLIP  for  May  9-10 ,  1980.  Predicted 
profiles  are  for  the  Garwood ,  Price ,  and  TL  models.  The  horizontal  scale  is  1  day  —  1.43  °C. 
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